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ABSTRACT 

A quadratic nonlinear generalization of the linear Rotta model for the slow pressure- 
strain correlation of turbulence is developed. The model is shown to satisfy realizability 
and to give rise to no stable non-trivial equilibrium solutions for the anisotropy tensor in 
the case of vanishing mean velocity gradients. The absence of stable non-trivial equilibrium 
solutions is a necessary condition to ensure that the model predicts a return to isotropy for 
all relaxational turbulent flows. Both the phase space dynamics and the temporal behavior 
of the model are examined and compared against experimental data for the return to 
isotropy problem. It is demonstrated that the quadratic model successfully captures the 
experimental trends which clearly exhibit nonlinear behavior. Direct comparisons are also 
made with the predictions of the Rotta model and the Lumley model. 
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1. Introduction 


The return to isotropy problem in turbulence has served as a cornerstone for the calibration 
of models for the slow pressure-strain correlation. The slow pressure-strain correlation is 
the part of the pressure-strain correlation that is independent of the mean velocity gradi- 
ents. It is experimentally observed that an initially anisotropic, homogeneous turbulence 
produced by the application of constant mean velocity gradients undergoes a relaxation 
to a state of isotropy when the mean velocity gradients are removed. Rotta 1 was the first 
to develop a turbulence model for the slow pressure-strain correlation which captured the 
return to isotropy behavior within the framework of second-order closure models. While 
the model (which is linear in the anisotropy tensor) performs reasonably well for small 
initial anisotropies, it can give rise to considerable errors for more general turbulent flows 
undergoing a relaxation from an initial state that is strongly anisotropic. This was pointed 
out by Lumley 2 who emphasized the need to incorporate nonlinear effects and developed a 
general nonlinear representation for the slow pressure-strain correlation based on invariant 
tensor representation theory. Subsequent experiments (see Gence and Mathieu 3 and Choi 
and Lumley 4 ) have unequivocally confirmed the nonlinear nature of the return to isotropy 
problem. 

In this paper, a quadratic model for the slow pressure-strain correlation, which contains 
only one independent constant, will be derived and calibrated. It will be demonstrated 
that this simple nonlinear model provides for a more complete description of the return to 
isotropy problem than either the linear Rotta model 1 or the quasi-linear model of Lumley 2 
and Shih and Lumley 5 (where nonlinearity is introduced through the invariants of the 
anisotropy tensor). In addition to doing a fairly good job in matching the results of 
four independent sets of experimental data, the quadratic model will be shown to satisfy 
realizability and to possess no stable non-trivial fixed points for the anisotropy tensor in 
relaxational turbulent flows. The latter constraint, which appears to have been overlooked 
in the previous literature on the subject, is a crucial constraint that nonlinear models 
must satisfy in order to ensure the return to isotropy in relaxational turbulent flows. 
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The mathematical properties and dynamical performance of the proposed model will be 
discussed in detail. 

2. Formulation of the problem 

In conventional, second-order turbulence modeling, transport equations for the Reynolds 
stress tensor u?uj and the dissipation rate e are solved in addition to solving the Reynolds 
averaged Navier-Stokes equations for the mean velocity U{. The exact transport equation 
for the Reynolds stress is 


dfUjUj + U k (uiUj) ik — Pij Tij kik Dij -f- n ti 


( 1 ) 


where 


Pij — UiU k U j ik UjU k U i ik 

T ijk = UjUjUjt + ~[pu;6j k + pujhk] - i/(u^7) ifc 

P 

Dij — 

% = ^iWU+Wjl) ( 2 ) 

given that u,- is the fluctuating velocity, p is the fluctuating pressure, p is the mass density, 
and v is the kinematic viscosity. In (1), P tJ is the production, Tij k is the diffusive trans- 
port, is the dissipation rate tensor, II tJ is the pressure-strain correlation, an overbar 
represents an ensemble mean, and a comma denotes a partial derivative with respect to 
the spatial coordinates. A transport equation is also usually carried for the dissipation 
rate e, which is defined as follows 

e = i/UijUij (3) 

The typical approach to modeling the unknown correlations T ijk , and II tJ in (1) is to 
express them as functions of the Reynolds stress tensor uiuj, the gradients of the Reynolds 
stress tensor (u.u/),*, the mean velocity gradients Uij, the dissipation rate e, and the 
kinematic viscosity v at lower turbulence Reynolds numbers. 
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At high Reynolds numbers, the dissipation rate tensor is nearly isotropic so that we 
can set 

A, = (4) 

The quantity IT,,- can be decomposed into a slow pressure-strain term II?- which is inde- 
pendent of Uij and a rapid term 11^ which is linear in Uij . The dependence of the slow 
pressure-strain correlation IT^- on u,it ; and e can be studied in isolation by considering 
the return to isotropy problem. The return to isotropy problem refers to the flow wherein 
homogeneous, anisotropic turbulence produced by mean velocity gradients is observed to 
relax towards isotropy when the mean velocity gradients are removed. This return to 
isotropy problem is an idealized case, but it is important nevertheless because it serves to 
extract and calibrate a part of the general functional dependence of n tJ - for more complex 
turbulent flows. 

Eq. (1) for the Reynolds stress tensor simplifies considerably for the return to isotropy of 
homogeneous turbulence since Uij and T, are identically zero, and Z),y may be assumed 
to be isotropic at high turbulence Reynolds numbers. The simplified form of (1) is written 
as follows: 

W7= ng - !*<« (s) 

where we have made use of the fact that the rapid part of the pressure-strain correlation 
vanishes for vanishing Uij. 

The anisotropy tensor is defined by 


bii = 


UiUj 






( 6 ) 


where 

q 2 = u,u< 


(7) 


The matrix associated with the anisotropy tensor is denoted by b. Using (5) and its 
contraction, the following evolution equation for bj is obtained 


bit = + 2-Ui 


9 2 ? 2 


( 8 ) 
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By defining a dimensionless tensor $,-y by 


rig. = 


( 9 ) 


(8) is rewritten as follows, 

bij — — — 2&tj) (10) 

Eq. (10) must be solved along with evolution equations for g 2 and c. For the return to 
isotropy problem, the exact transport equation for q 2 and the standard model equation for 
e are as follows 

q 2 = —2e 

£ = -2 C (2 ~ (11) 

Q 2 

where C f2 is a constant that assumes a value of approximately 1.90. 

The important variables in this homogeneous, high Reynolds number flow are 6,y, q 2 
and e. Lumley 2 gave the general tensorial form of the function $,y(6,y, q 2 , e) which can be 
written as follows 

= ai(//, III)bij + a 2 (II,III){b ik b kj - y%] (12) 

In (12), c*i and a 2 are general functions of the invariants II = b ik b ki and III — b ik b k jbji. 
It should be noted that II and III denote the traces of b 2 and b 3 respectively and differ 
from the principal invariants of b by constant multiplicative factors since b is traceless. 
On substituting the model of $,y from (12) into (10), we obtain 

bij = - 4 { [«!(//,///) - 2 It,,- + «,(//, III)[bi t b ki - I } (13) 

q o 

A prescription of the functions aj (//,///) and a 2 (II, III) is now required to close 
(13). The approach that is most widely used in practice is the linear model proposed by 

Rotta 1 in which cti (II, III) is assumed to be a constant and a 2 (II,III) is assumed to 

be zero. The commonly used second-order closure model of Launder, Reece and Rodi 6 is 
based on the Rotta model where 


«!(//,///) = 0 ^ 3.0 

a 2 (II,III) = 0 (14) 


4 


Lumley and Newman 7 and Lumley 2 have pointed out that experimental data does not 
support the linear model of Rotta. The problem with the linear model is twofold: (a) its 
prediction that each component of the anisotropy tensor decays at the same rate, and (b) 
its prediction that this decay rate is independent of the initial state of anisotropy. These 
predictions can be seen mathematically in the closed form solution of the return to isotropy 
problem for the Rotta model, given by 

6„ = MO) |1 + 2(0., - (15) 

9o 

In contradiction to (15), experiments indicate that the rate of the return to isotropy 
decreases with increasing values of the third invariant III, and that different components 
of the anisotropy tensor relax to zero at different rates. These are effects that are decidedly 
nonlinear. 

In order to incorporate nonlinear effects, Lumley 2 developed a model within the frame- 
work of (12), and also allowed the coefficients ai and a 2 to additionally depend on the 
turbulence Reynolds number Ret = q*/9eu. This quasi-linear model of Lumley 2 is as 
follows 

a^IIJII.Re,) = 2 .0+|exp(^){-^= 

+80.1 ln[l + 62.4 (—IT + 2.3///*)] } 

a 2 {II,III,Re t ) = 0 (16) 

where 

F = l + 9/r + 27//r 

//' = -y, III' = ^ (17) 

Here, we refer to this model as being quasi-linear since there is no tensorial nonlinearity, 
that is, 

$ij — oci(II,III,Re t )bij (18) 

It will be shown later that while this model does correctly mimic some of the nonlinear 
effects associated with the return to isotropy (see Shih and Lumley s ) it yields an incorrect 
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representation of the problem in phase space. In the next section, a quadratic model will 
be developed that appears to capture the essence of the nonlinearity in phase space. 

3. The quadratic nonlinear model 

We will now develop a quadratic model, which is a special case of (12). Since each com- 
ponent energy (where Greek indices indicate no summation) is bounded by 

0 < u a u a < q 2 (19) 

it is a simple matter to show from (6) that each principal value of 6,-y (denoted by ) 
lies in the range 

- | < b {a) < l ( 20 ) 

Consequently, the principal values of are less than unity, and it is possible that a 
lower order Taylor expansion of $,y(6,y) around the isotropic state of 6,y = 0 can yield an 
adequate approximation. Carrying out such a Taylor expansion around 6,y = 0 , while 
constraining to be of zero trace, gives 

* = C x b-C 2 (b 2 -yI) + C 3 (b 3 -^I) + ... (21) 

where C x , C 2 , C 3 , ... are constants. In (21) we use coordinate-free notation, that is, we 
denote 6«y and 6 { j by the associated matrices <£, b and I. Using the Cayley-Hamilton 
theorem to express b 3 and higher powers of b in terms of II, III, b and b 2 reduces 
(21) to (12) where oq [II, III) and a 2 (II, III) are polynomial functions of II and III of 
infinite degree. 

The model of that we propose in this paper is a quadratic truncation of (21), and 
is given by 

= Cxbij - C 2 [b ik b kj - y*,-] (22) 

where C\ and C 2 are constants. Such a quadratic model is attractive because it is the 
simplest nonlinear model for (of course a linear truncation of (21), which leads to 
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Rotta’s model, would be even more simple but, unfortunately, is in conflict with experi- 
mental data) . Furthermore, since for most problems of engineering interest the magnitude 
of bij tends to be less than 1/3, there is a good chance for a quadratic model to yield an 
acceptable approximation for a variety of flows. 

At this point, we will derive the evolution equations determining the phase space dy- 
namics for the return to isotropy problem. It is useful to introduce the transformed time 
variable r defined by 

dr = ~dt (23) 

<r 

in order to free the formulation of the problem from any explicit dependence on c and 
q 2 . Using the system of equations (11) to solve for q 2 (t) and e(t), and integrating (23) we 
obtain the following explicit relationship between r and t : 

+ < 2<) 

Using (23), the transport equation (13) for 6,-y can then be rewritten as 

^ = -{Mil, III) - 2 )bij + o 2 (II, III) [bi , k b ki - } (25) 

which is clearly independent of q 2 and e. 

We now show that the phase space for ( 25 ) is two-dimensional. Let the principal axes 
of b^ at r = 0 be {xi,x 2 , x 3 } and the corresponding eigenvalues of 6,,- be denoted by b n , 
b 22 , and 6 33 . The off-diagonal terms, (t ± j) are zero at r = 0. On examining (25) in 
principal axes form it is clear that 

- — = 0 for i 7^ j (26) 

dr 

Since 6»y(0) is zero for t 7^ j, we conclude from (26) that the off-diagonal terms of 6 tJ are 
identically zero for all times; and therefore, the principal axes of do not change in time. 
The time invariance of the principal axes has been experimentally observed, for example 
by Gence and Mathieu 3 . Therefore the only time varying components of 6,y, referred to the 
coordinate system {xi,x 2 ,xz}, are bn, 6 22 , and 633. Since 6 tJ is traceless, only two of the 
three quantities 6n, 6 22 > and 633, are independent. Thus the phase space in this problem is 
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two-dimensional. The two phase space variables may be chosen to be any two independent 
functions of the eigenvalues of 6,y. 

A possible choice for phase space variables is II and III. The variables II and III 
are desirable to use since they are uniquely determined by the principal values of are 
invariant to a rotation of coordinates, and vanish in isotropic turbulence. Furthermore, 
because II > 0 the phase space becomes semi-infinite in the variable II. The evolution 
equations for J/(r) and III(t) will now be derived. Multiplying (25) by 6,y and babkj, 
respectively, we obtain the following equations 

^=-2{(a, -2 )II + a,III} (27) 

^ = — 3{ (a, - 2)/// + a=a[tr (b 4 ) - ~\ } (28) 

In (28),tr (b 4 ) denotes the trace of the tensor b 4 . Using the Cayley-Hamilton theorem, it 
follows that 

b 4 = I//b ! + ™b (29) 

By taking the trace of (29) , and substituting the result into (28) , we obtain 

^ = — 3{ (a x - 2)1 1 1 + a 2 ~} (30) 

Equations (27) and (30) constitute a coupled nonlinear system of ordinary differential 
equations for 7/(r) and IJJ(r). 

Both the Rotta model (14) and the Lumley model (16) choose a 2 = 0. When a 2 = 0, 
(27) and (30) can be combined and integrated to yield the general solution 

(II) 1 ' 2 = c(JII) 1 / 3 (31) 

where c is a constant of integration determined by initial conditions. In terms of the 
variables 


i = (my/ 3 

v = (II) 112 (32) 
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( 33 ) 


which were first introduced by Choi and Lumley 4 , (31) takes the linear form 

Vo £o 

Thus, these alternative variables have the advantage that in the £ — rj phase space both the 
linear and quasi-linear models show straight line trajectories. Some of the experimental 
data that will be shown later exhibit curved trajectories in the £ — rj phase space which 
constitutes conclusive evidence in support of the need for nonlinear models. 


4. Fixed point analysis in phase space 

A fixed point analysis of the quadratic model in the phase space of the invariants of 6,-y 
will now be presented. The fixed points of (27) and (30) (where, for the quadratic model, 
c*i = Ci and a 2 = ~&2 ) are the equilibrium solutions in the limit as r — ► oo which are 
denoted by (//<„, 1 1 loo)- If the turbulence is to return to isotropy 

I loo = ///oo = 0 (34) 


Therefore, since physical and numerical experiments indicate that any initial anisotropy 
which is realizable tends to relax to zero when the mean velocity gradients are removed, 
it is essential that nonlinear models for the slow pressure-strain correlation have no stable 
fixed points other than (0,0). 

Since the linear model (for which C 2 = 0) gives rise to the following system of evolution 
equations 


dH_ 

dr 

dm 

dr 


- 2 {Ci - 2)11 
-3 (Ci - 2)/// 


(35) 


it is a simple matter to show that (0, 0) is the only fixed point, and that it is a stable node 
provided that Ci > 2. 

We now examine the existence and stability of the fixed points predicted by our 
quadratic model (22). The quadratic model leads to the following equations for //(r) 
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and III(t): 


— = —2(Ci — 2)11 + 2C 2 III 
cLt 

^ (36) 

The fixed points of the system (36) are solutions of the nonlinear algebraic equations 
obtained by setting dll/dr = dill /dr = 0, and are given by 

7/00 = 0 , 1 1 loo = o (37) 

Hoo = 6/3 2 , II loo = 6/3 3 (38) 

where /? = [C\ — 2)jC 2 . In order to examine the stability of these fixed points we introduce 
the new variables (zi,z 2 ) defined by 

Zl = II- I loo , *2 = III - 1 1 loo (39) 

If we substitute (39) into (36), and retain only the terms that are linear in z\ and z 2 , 
we obtain the following linear system of ordinary differential equations 

dzi . , v 

— - AijZj (40) 

where is the following matrix 

-2(C' 1 - 2) 2 C 2 

C 2 II 00 -3 (C, - 2) 

The eigenvalues of Aij determine the stability of the fixed points. A fixed point is stable if 
all the eigenvalues of A, ; are negative. The fixed point (0, 0) has eigenvalues A^ given by 

A^ = —2{Ci - 2), A (2) = —3 [Cx - 2) (41) 

and, hence, is a stable node provided that 

C x > 2 (42) 
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On the other hand, the fixed point (6/? 2 ,6/? 3 ) has the eigenvalues 

AW = Ci- 2, A< 2 > = -6(0! - 2) (43) 

Since the eigenvalues are of opposite sign no matter what value is taken for Ci, the fixed 
point (6/? 2 , 6/? 3 ) is a saddle which is unstable . It is thus clear that IIIoo — I loo — 0 is the 
only stable equilibrium solution and the quadratic model predicts a return to isotropy for 
all initial conditions on 6,-y. 

5. Constraints on the model constants 

The quadratic model has two constants Ci and C 2 which, for the purpose of stability, 
need only satisfy the constraint Ci > 2. We will now explore whether there are additional 
physical constraints on the range of allowable values of C 1 and C 2 . Choi and Lumley 4 and 
Le Penven, Gence and Comte-Bellot 8 have analyzed available sets of experimental data 
(approximately ten in number) on the relaxation of anisotropic turbulence. They conclude 
from the data that a characteristic nondimensional time rate for the relaxation (defined 
by decreases with increasing values of III. An examination of (36) shows that in 

order for the model to reproduce this trend we must choose 

C 2 > 0 (44) 

Lumley 2 uses the concept of realizability, which requires that model predictions satisfy 
all moment inequalities; for example, the model is required to yield non-negative compo- 
nent energies. Lumley’s method leads to various equations relating the model constants. 
Pope 9 suggests that a related but weaker mathematical statement, which constitutes a 
sufficient condition to guarantee non-negative component energies, should be used. The 
application of both approaches in the context of the quadratic model for <&,•/ will now be 
discussed. 

The transport equations for 6 fJ - and u,u ; in the present problem are 

= 26«) (45) 
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The three eigenvalues of ujuj are the three components in the principal axes of ujuj ; 
these components are u 2 , u\ and uf. Consider one of the eigenvalues, for example u\. 
Realizability requires that 

«! > 0 (47) 

Using the Lumley 2 approach to satisfy realizability leads to the requirement that 

u\— 0 when u\ = 0 (48) 

whereas the Pope 9 approach leads to the requirement that 

u\> 0 when u 2 = 0 (49) 


Both approaches guarantee that the realizability condition (47) is satisfied. However (49) 
can prevent uf from ever becoming zero (which is the physically realizable limit of two- 
dimensional turbulence), and hence is more conservative than (48) in this regard. But 
since (49) does guarantee positive component energies, we will be content for the moment 
to consider it since it places less restrictions on C t and C 2 . 

Since (49) is to be satisfied, using (46) we require that 


— e 4>n — -e > 0 when u\ = 0. 

O 


(50) 


This implies that 


4> u < — - when u? = 0 
3 


Substituting the quadratic model for from (22) into (51) leads to the inequality 

Cib n - C 2 (b 2 n - y) < when u\ = 0 
Since b n = —1/3 when u\ = 0, (52) becomes 


(51) 


(52) 


Ci - 2 > C 2 {II - -) 


(53) 


From the constraint (42) , it is clear that the left-hand-side of (53) is positive. Since from 
(44) C 2 > 0, the maximum value of the right-hand-side of (53) occurs when II is at its 



maximum value. Therefore a sufficient condition for satisfying (53) is obtained when II is 
equal to its maximum value II m ax • It can be easily shown that 

IU, = \ (54) 

By substituting the maximum value of JJ, which is given by (54), into (53), we obtain the 
following constraint 

C 2 < 3 (Ci - 2) (55) 

For a given value of Ci, (55) determines the maximum allowable value of C 2 . 

6. Model performance 

The coupled system of nonlinear ordinary differential equations (27) and (30) were nu- 
merically integrated by a Runge-Kutta method to obtain II{r) and III(r). The initial 
conditions corresponding to those of the previously conducted experiments were used. The 
phase plots in £ — rf space and the transient plots of II(t) were the yardsticks by which 
the model results and experimental data were compared. 

A numerical optimization showed that for a given value of C\ it was best to choose the 
maximum value of C 2 allowed by (55), that is, to set 

C 2 = 3(Ci - 2). (56) 

Due to (56), C 2 is uniquely determined by C\ making the model one with a single free 
constant C\. It is interesting to note that (56) is similar to the constraint that is obtained 
by using Lumley’s method of enforcing realizability. Thus it appears that Lumley’s method 
works better for this problem. The following choice of the one independent model constant 

Ci = 3.4 (57) 

seems to best reflect the experimental data. 

Results on the performance of the linear Rotta model, quasi-linear Lumley model, and 
the present quadratic model are compared with four sets of experimental data in Figs. 1-8. 
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Figs. 1-2 consider the plane distortion experiment of Choi and Lumley 4 . The experimental 
data in Fig. 1 shows strong non-linearity in the £ — rj phase plane (as demonstrated by a 
curved trajectory), which is reproduced by the quadratic model. However, the trajectories 
in £ — r/ phase space of the Rotta and Lumley models are always straight lines, and are 
therefore in serious disagreement with this experimental data. Fig. 2 shows the temporal 
decay of II. It is seen that the present quadratic model is superior to both the Rotta model 
and the Lumley model with regard to capturing the temporal behavior of II. Figs. 3-4 
pertain to the axisymmetric expansion experiment of Choi and Lumley 4 . The experimental 
phase space trajectory seems to be fairly linear as shown in Fig. 3, and the present quadratic 
model also shows almost linear behavior. Apparently, the quadratic model not only shows 
nonlinear phase trajectories in £ — rj phase space when called for, but also gives rise to 
approximately linear trajectories in the appropriate region of phase space where nonlinear 
effects are small. It is seen in Fig. 4 that the results of both the Rotta model and the 
quadratic model on the temporal behavior of II are in reasonable agreement with the 
experimental data; however, the Lumley model predicts a decay rate of II which is much 
more rapid than the experimentally observed result. 

Le Penven, et al. 8 used two different plane contractions to generate anisotropic turbu- 
lence with positive III in one case and negative III in the other case. Their experimental 
data on the relaxation of turbulence with positive III is shown in Figs. 5-6. Their data of 
Fig. 5 indicates a moderately curved trajectory in phase space, and the quadratic model 
yields a prediction that is in qualitatively good agreement with the data. Fig. 6 shows that 
the quadratic model does a better job than the Lumley model in predicting the temporal 
decay of II, and that the Rotta model is slightly better than the quadratic model with re- 
gard to predicting the temporal behavior of II. In Figs. 7-8, we consider the experimental 
data of Le Penven, et al. 8 pertaining to the case with negative III. The phase portrait 
of Fig. 7 shows that their experimental data exhibits mild curvature in the £ — ri phase 
space, which is reproduced by the quadratic model. Fig. 8, which depicts the temporal 
evolution of II, shows that the Rotta model grossly underpredicts the rate of decay of II. 
The Lumley model is seen to be somewhat better than the proposed quadratic model in 
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capturing the rapid decay of II in this experiment. Nevertheless, it is clear from the results 
presented in this section that the quadratic model does the best overall job in predicting 
the experimental trends. 

7. Conclusion 

A quadratic generalization of the linear Rotta model has been developed for the slow 
pressure-strain correlation of turbulence. This simple nonlinear model, which has only one 
independent constant, was shown to constitute an improvement on both the Rotta model 
and the quasi-linear Lumley model in predicting the trends observed in return to isotropy 
experiments. Furthermore, the quadratic model was shown to satisfy realizability and to 
have no stable non-trivial equilibrium solutions for the anisotropy tensor in relaxational 
flows- a crucial constraint which guarantees that the model predicts a return to isotropy, 
independent of the initial conditions. The most notable feature of this new quadratic model 
is its ability to capture the essence of the return to isotropy problem with a nonlinear 
structure that is substantially simpler than previously proposed nonlinear generalizations 
of the Rotta model. 

Future research can be directed on two fronts: the development of a low turbulence 
Reynolds number correction of the model and the addition of successively higher order (for 
example, cubic) nonlinear terms. Insofar as the former research direction is concerned, the 
direct numerical simulations of Lee and Reynolds 10 have indicated that history dependent 
effects of the Reynolds stress and dissipation rate anisotropy tensors are important in low 
Reynolds number turbulent flows. These are complicating features that are not likely to be 
described by simple nonlinear models with an algebraic structure. In regard to the addition 
of higher degree nonlinearities, we feel that it is best to first provide a systematic test of 
a quadratic model since second-order closures that are quadratic in the anisotropy tensor 
can satisfy constraints such as material frame-indifference in the limit of two dimensional 
turbulence 11,12 and realizability. The satisfaction of these constraints, which are violated 
by simple second-order closures such as the Launder, Reece and Rodi model 6 , allow for the 
description of more complex turbulent flows within the framework of a workable theory. A 
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systematic program for the development of improved second-order closure models that are 
quadratically nonlinear is currently underway and will be the subject of a future paper. 
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